Le réplicat WT.K1 est différent des autres : il présente une forte hétérogénéité, avec soit de très forts ou très faibles counts suivant les gènes. Après normalisation, il semble toujours se démarquer des autres et fausser l’analyse d’expression différentielle.
Ce document résume les différences entre les design avec ou sans ce réplicat, et les résultats obtenus avec DESeq2 et EdgeR.
Aspect du réplicat moche
load("y_all.RData")
yAll = y
load("y_minK1.RData")
yMinK1 = y
normAll <- cpm(yAll, normalized.lib.sizes = TRUE)
normMinK1 <- cpm(yMinK1, normalized.lib.sizes = TRUE)
not_norm <- cpm(yAll, normalized.lib.sizes = FALSE)
rd_genes = sample(rownames(normAll), 500)
heatmap(not_norm[rd_genes, ], main = "cpm Sans normalisation")
On voit que la matrice normalisée est beaucoup plus homogène sans K1. Même conclusions avec la méthode de normalisation de DESeq2 non montrée ici.
Comparaison des log Fold Change
Ici, avec EdgeR, on va comparer les logFoldChange entre les deux designs.
load("et_all.RData")
etAll = et
load("et_minK1.RData")
etMinK1 = et
plotMD(etAll)
abline(h = c(-1, 1), col = "blue")
abline(h = c(0), col = "red")
2-1
Down 27
NotSig 5832
Up 75
2-1
Down 78
NotSig 6573
Up 151Le log fold change n’est pas centré en 0, probablement car le nombre de counts de K1 est plus fort que les autres, et tire le ratio de manière systématique vers une valeur inférieure à 1.
Comparaison des gènes détectés
load("DEGenesEdgeRAll.RData")
EdgeAll = DEGenesNames
load("DEGenesEdgeRWithoutK1.RData")
EdgeminK1 = DEGenesNames
load("DEGenesDESeqAll.RData")
deseqAll = DEGenesNames
load("DEGenesDESeqWithoutK1.RData")
deseqminK1 = DEGenesNamesOn détecté déjà beaucoup plus de gènes avec EdgeR, différences dues à leur méthode de normalisation et probablement à leur estimation de la variance et tests de significativité.
df <- data.frame(matrix(ncol = 2, nrow = 4))
x <- c("Design", "GeneNumber")
colnames(df) <- x
df$Design <- c("DESeq", "DESeqMinK1", "EdgeR", "EdgeRMinK1")
df$GeneNumber <- c(length(deseqAll), length(deseqminK1), length(EdgeAll), length(EdgeminK1))
ggplot(data = df, aes(x = Design, y = GeneNumber)) + geom_bar(stat = "identity",
fill = "steelblue") + ggtitle("Nombre de gènes différentiellement exprimés suivant le design et le package R") +
geom_text(aes(label = GeneNumber), vjust = 1.6, color = "white", size = 3.5) +
theme_minimal()
Les tests permettent de détecter plus de gènes lorsque le réplicat douteux est retiré de l’étude.
x <- list(A = deseqAll, B = deseqminK1, C = EdgeAll, D = EdgeminK1)
ggVennDiagram(x, category.names = c("DESeq", "DESeqMinK1", "EdgeR", "EdgeRMinK1"))
[1] "NRT2.1" "GA3OX1" "TCP23" "NIA2" "SDC" "PGD1" "HHO2"
[8] "NIA1" "BXL7" "NIR1" "AHK1" "FD3" "PXC1" "MDH"
[15] "BT2" "LBD38" "BAT1" "NLP3" "PUB28" "G6PD2" "AO"
[22] "BT1" "LBD37"